Equilibrium Properties of Double-Screened-Dipole-Barrier SINIS Josephson Junctions 



Branislav K. Nikolic, 1 J. K. Freericks, 1 and P. Miller 2 
1 Department of Physics, Georgetown University, Washington, DC 20057-0995 
2 Department of Physics, Brandeis University, Waltham, MA 02454 

We report on a self-consistent microscopic study of the DC Josephson effect in SINIS junctions 
where screened dipole layers at the SN interfaces generate a double-barrier multilayered SIN struc- 
ture. Our approach starts from a microscopic Hamiltonian denned on a simple cubic lattice, with 
an attractive Hubbard term accounting for the short coherence length superconducting order in the 
semi-infinite leads, and a spatially extended charge distribution (screened dipole layer) induced by 
the difference in Fermi energies of the superconductor 5* and the clean normal metal interlayer N. 
By employing the temperature Green function technique, in a continued fraction representation, the 
influence of such spatially inhomogeneous barriers on the proximity effect, current-phase relation, 
critical supercurrent and normal state junction resistance, is investigated for different normal in- 
terlayer thicknesses and barrier heights. These results are of relevance for high-T c grain boundary 
junctions, and also reveal one of the mechanisms that can lead to low critical currents of apparently 
ballistic SNS junctions while increasing its normal state resistance in a much weaker fashion. When 
the N region is a doped semiconductor, we find a substantial change in the dipole layer (generated 
by a small Fermi level mismatch) upon crossing the superconducting critical temperature, which 
is a new signature of proximity effect and might be related to recent Raman studies in Nb/InAs 
bilayers. 

PACS numbers: 71.27.+a, 74.50.+r, 74.80.Fp, 73.40.Jn 



I. INTRODUCTION 

The Josephson effectHl is one of the most spectacular 
phenomena arising from the macroscopic phase coher- 
ence of Cooper pairs. A dissipationless current flows at 
zero voltage between two superconductors weakly cou- 
pled through a tunnel barrier (SIS, where S and I de- 
note a superconductor and an insulating barrier, respec- 
tively) or weak links (ScS, SNS, etc., where c stands 
for a constriction, and N for a normal metal). The 
study of such inhomogeneous superconducting structures 
has been driven by both interest in the fundamentals of 
quantum mechanics, and by the potential application of 
JosepJison junctions as circuit elements in electronic de- 
vices. □ 

Recently, considerable attention hasJpeen directed to- 
ward the study of SINIS junctionsJ3~Q where the insu- 
lating tunnel barrier is split into two pieces separated by 
a normal metal. These types of junctions have provided 
a playground to study the interplays between the meso- 
scopic coherence of a single-particle wave function in the 
normal metal and the macroscopic coherence of a many- 
body wave function of Cooper pairs.EI Furthermore, the 
reexamination of various multilayered structures of the 
SINIS type in applied research has been driven by the 
struggle to optimize the performance of Josephson junc- 
tions in loWj-temperature superconducting (LTS) digital 
electronics JErEJ In mesoscopic superconductivity, one fre- 
quently deals with S-Sm-S junctionsB (Sm being a heav- 
ily doped semiconductor with a two-dimensional electron 
gas) where the role of the I layer is played by a space- 
charge layer arising at the S-Sm interface (additional 
scattering at the interface can occur from the mismatch 



between the effective electron masses and Fermi momenta 
in the S and Sm). The technological advances in fabri- 
cating such hybrid structures^ have given an. impetus to 
the field of mesoscopic superconductivityuu where the 
two-dimensional electron gas is amenable to an engineer- 
ing of its "metallic" properties; i.e, one can tune the 
Fermi wavelength, or mean free path, and one can confine 
electrons with gate electrodes. In such structures, phasp. 
coherence of the electron and Andreev-reflected holctil 
at the SN interface can be studied without too much 
normal reflection, because the charge-accumulation layer 
arising at a typical Nb/InAs interface, or the Schottky 
barrier at a Nb/Si interface, are much-more transparent 
than typical dielectric tunnel barriers.u 

While initial understanding of the Josephson effect 
came from studies of tunnel junctions,!!! further devel- 
opments concentrated on weak linkaL3 which provide the 
non-hysteretic (i.e., single valued) I -^V characteristic 
needed for applications, like SQUIDsE^I or rapid single 
flux quantum logic.H-3 The return to SIS junctions came 
after the fabrication of Nb/Al tunnel junctionals with a 
reliable control of the critical current (conventional tun- 
nel junctions can be made non-hysteretic by externally 
shunting their high capacitance— with a resistor, which 
reduces the overall performancaia) . The renewed inter- 
ests in SINIS multilayered junctions for LTS electron- 
ics comes from an attempt to combine the advantageous 
properties of both weak links and tunnel junctionsu — the 
SINIS junctions are intrinsically shunted, while exhibit- 
ing large characteristic voltages with moderate critical 
current densities (in fact, rapid single flux quantum de- 
vices requite large critical current densities, to reduce the 
error rate,Efl which is difficult to achieve using standard 
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Nb/Al/AlO^Nb tunnel junction technology, but might 
be reachednin SINIS junctions with carefully engineered 
properties!^) . When the N interlayer is clean, the junc- 
tion resistance is mainly controlled by scattering at the 
interfaces (like in conventional Nb/Al/A10 x /Al/Nb junc- 
tionsO), and not by the interlayer material properties. 

Here we undertake a study of a special class of SINIS 
junctions where the double-barrier structure arises from 
two inhomogeneous screened dipole layers (SDL) deter- 
mined by a relatively large Debye screening length Id 
of a few lattice spacings. We start from a microscopic 
lattice Hamiltonian with the S and N layers described 
by different metals that have the same bandwidth, but 
their Fermi levels are misaligned. The Fermi level mis- 
match forces a charge redistribution, with the strongest 
deviation from uniformity located near the SN inter- 
face, which is gradually diminished inside the bulk lay- 
ers on a length scale set by Id- The charge profile en- 
sures an equilibration of the chemical potential through- 
out the system when no bias voltage is applied. Since 
we assume a screening length of a few lattice spacings, 
the dipole layer is spatially extended (i.e., thicker than 
just one monoatomic layer). This choice of microscopic 
junction parameters allows us to examine the charge 
redistribution appearing between conductors which are 
less efficient in screening than ordinary metals (such as 
the underdoped cuprates or InAs). Our treatment of 
the double SDL barrier is fully microscopic and self- 
consistent, meaning that effects of the static electric po- 
tential (generated by the excess charge) on the Joseph- 
son current and on the normal state resistance are re- 
lated to the parameters of the underlying Hamiltonian, 
rather than ^hapacterizing the barrier by an effective 
transparencyD'HEJ D, or using a delta function poten- 
tial |-a,t— the SN interface to model the normal reflec- 
tionc2lE3 (in addition to the inevitable retroreflcctionM) . 
We tackle both the fundamental aspects of the problem 
(like the self-consistent evaluation of the order param- 
eter, the change of its phase across the junction, and 
the emergence of non-sinusoidal current-phase relations) 
and issues relevant for applications (like the characteris- 
tic voltage, a product of the critical current I c and the 
normal state resistance Rn, which determines the high- 
frequency performance of the junction). Our junctions 
are three-dimensional (3D) and clean, so that quasipar- 
ticle transport through the N interlayer is ballistic. 

Previous theoretical work on ballistic SI NIS junctions 
focused on— .resonant supercurrents in low-dimensional 
structures£3~Ea Mesoscopic superconductivity coherence 
effects in 3D junctions (e.g., a current proportional to 
D of the barrier, rather than the characteristic D 2 de- 
pendence for two uncorrelated sequential tunneling pro- 
cesses) have been investigated in Ref. |]. These junctions 
are mostly similar to the ones studied here, except that 
our "microscopic" charge accumulation barriers are not 
atomically sharp interfaces that can be described by a 
phenomenological transparency D. A more microscopic 
treatment of the effect of charge inhomogeneity for 




FIG. 1. Microscopic stacked planar geometry of a Joseph- 
son junction defined on an infinite simple cubic lattice with a 
lattice constant a. The normal interlayer contains Nm planes 
(ranging from 1 to 60) which are coupled to semi-infinite su- 
perconducting leads (the junction thickness is L = Nno). 
These layers, together with the first Ns planes (30 in our cal- 
culations) in each lead, comprise the region of the junction 
where the self-consistent calculation is performed. The junc- 
tion is allowed to have spatial inhomogeneity only within the 
2Ns + Nm modeled planes, but the calculations are for an in- 
finite system. The insulating barriers are formed by a charge 
redistribution that is localized near the SN interfaces. 

normal transport through the contact of two differ- 
ent metals (a problem frequently appearineJn the mul- 
tilayers of giant magnetoresistance devicesus) has been 
undertaken using the Boltzmann equationpD and in su- 
perconducting junctions using-.quasiclassical methods in 
a non-self-consistent fashion.Ea It is worth emphasizing 
that standard quasiclassical Green function techniques, 
which exploit the fact that macroscopic quantities vary 
on a length scale substantially exceeding the interatomic 
distance, cannot be applied directly to problems contain- 
ing boundaries between two different metals. Since elec- 
tron reflections lead to fast spatial variations of the orig- 
inal Green functions around the boundary, the method 
has to be extended properly to take this into account (see 
Ref. |28| for details) . 

Our study is relevant for three types of recently ex- 
plored, experimental systems: (i) grain boundary junc- 
tionsEj in high-T c superconductors, where our short co- 
herence length superconductor and the poor screening 
of the excess charge (i.e., Debye screening length com- 
parable to the coherence length), mimic the effect of a 
charge imbalance at the grain boundaries on the depres- 
sion of the order parameter, and thereby the intergranu- 
lar current densityEj'Eil (without complicating the prob- 
lem further with d- wave symmetry); (ii) Raman scatter- 
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ing studies&2l of the proximity effects in Nb /InAs hybrid 
structures reveal a substantial change of the charge ac- 
cumulation layer formed at such interface above and be- 
low the T c of Nb — we also find that I layer induced by 
a small Fermi level mismatch is modified by proximity 
effects in our SINIS junctions when the carrier concen- 
tration in the N is 100 times smaller than in thc.5 1 ; (iii) 
recent experiments on ballistic SNS junctions,E3 in the 
limit where I c and Rn do not depend on the thickness of 
the N, exhibit a much smaller characteristic voltage than 
predicted for short clean ScS junctions — the scattering 
off a dipole charge layer is an example of a process which 
sharply reduces I c , but only weakly increases Rjy. 

The paper is organized as follows. In Sec. |l| we intro- 
duce the model and the main ideas of the Green function 
computational technique (employed to solve the quantum 
problem of the charge distribution and equilibrium trans- 
port; the electrostatic problem of the potential generated 
by these charges is solved classically). Section ID con- 
tains the results for the self-consistent pair amplitude (or 
the order parameter) and the local change of the phase 
across the junction. The current-phase relation for differ- 
ent strengths of the electrostatic potential generated by 
the SDL is discussed in Sec. IV, where we also evaluate 
the characteristic voltage I c Rn- We conclude in Sec. |v|. 



II. MODELING A SINIS JUNCTION WITH A 
DOUBLE-BARRIER SCREENED DIPOLE LAYER 

Early studies of the Josephson effect in SINIS junc- 
tions were based on a tunneling Hamiltonian formalisM. 
and perturbation theory in the barrier transmissiijity.tSI 
Later on, quasiclassical Green function techniquestJ were 
applied to a double-barxier junction with the N inter- 
layer in the dirty limit .El While these results are valid 
only in a few limiting cases, a recent reexamination of 
this problem covers a wider range of parameters.ula For 
example, when transport through the N interlayer is bal- 
listic (mean free path greater than ther-thickness of the 
junction), one cannot use standard toolset like the Usadel 
equation. Instead, a solution of the Gor'kov equations 
for the Green, functions of the double-barrier structure 
is required .lEj Furthermore, if the I barriers are not of 
low transparency, the usual-arguments for the validity 
of rigid boundary conditionstj (i.e., taking the gap A to 
be constant inside the superconducting leads) fail when 
the S and N regions have the same cross section, and 
the thickness of the junction is not much larger than the 
superconducting coherence length £ - ln such cases, the 
critical current density can be close to the bulk critical 
current density, and a self-consistent evaluation of the 
order parameter inside the S is needed to ensure cur- 
rent conservation throughout the structure.Ej~E!3 Since 
we choose to work with a short coherence length super- 
conductor, quasiclassical approximations neglecting dy- 
namics on a length scale below £o are not applicable (in 



our case £o ~ 4a is not much larger than the Fermi wave- 
length Xp w 2a, and spatial variation of the order pa- 
rameter A on a length scale smaller or comparable to £o 
is important). 

Our approach to quantum transport in ballistic 
SINIS junctions starts from a microscopic Hamiltonian. 
defined on a simple cubic lattice (of lattice constant a) .Eil 
It allows us to describe the transport for an arbitrary 
junction thickness, temperature, and barrier strength. 
Also, the geometry is such that the N interlayer has the 
same width as the S leads. For computational purposes, 
the infinite lattice which models the junction is divided 
into a self-consistent part and a bulk part, as shown in 
Fig. [l] A negative-t/ Hubbard term is employed to model 
the real-space pairing of elections due to a local instanta- 
neous attractive interaction.aO The lattice Hamiltonian 
is given by 



(1) 



where c\ a (c^) creates (destroys) an electron of spin a 
at site i, t^j is the hopping integral between nearest- 
neighbor sites i and j (energies are measured in units 
of t), which is taken to be the same in the S and N, 
and Ui < is the attractive Hubbard interaction for 
sites within the superconducting planes. The normal in- 
terlayer is described by the noninteracting part of the 
Hamiltonian (|l]), which is just a (clean) nearest-neighbor 
tight-binding model with a diagonal on-site potential V% . 
The potentials Vi are not given a priori, but instead are 
calculated self-consistently by first determining the local 
electronic charge density and comparing it to the bulk 
charge density of the corresponding S or N layers. The 
imbalanced charge on each plane generates an electric 
field and thereby an electric potential. Summing the 
contributions from the charges on all other planes then 
yields the total local potential Vf and the local potential 
energy shift Vi — eVf . We now recalculate the charge 
density on each plane and iterate until Vi is determined 
self-consistently (see below for a detailed description of 
the algorithm). The local potentials Vi are largest near 
the SN interface, and decay as one approaches the bulk 
leads. 

We use the Hartree-Fock approximation (HFA) for the 
interacting part of the Hamiltonian ([!]). This accounts 
for the superconductivity in the S region in a way which 
is completely equivalent to a conventional BCS theory 
with an energy cutoff determined by the electronic band- 
width rather than by the phonon frequency. We choose 
half-filling ns — 1 and Ui — —2 on the sites in the super- 
conducting leads. The homogeneous bulk superconduc- 
tor has a transition temperature T c = 0.11 and a zero- 
temperature order parameter A = 0.198. This yields a 
standard BCS gap ratio 2A/ (ksT c ) w 3.6 and a short 
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FIG. 2. Scaling of the local density of electrons with the 
thickness L — Nnci of SINIS Josephson junctions. The fill- 
ing is the same at each site within the planes (plane number 
one is the first plane along the z-axis inside the self-consistent 
region from Fig. |l|). The difference in Fermi energies of the 
S and N is AEf = 3.0. The bulk equilibrium value of 
the charge density is set to half-filling in both the S and N 
(ns = tin ~ 1)- The inhomogeneous charge redistribution is 
(approximately) symmetric around the SN interface only for 
a thick enough junction. In the case of a thin junction, the 
charge is depleted in the N interlayer to lie below half-filling 
since the screening length is a few lattice spacings. 

coherence length £o = hvp/(nA) ~ 4a. The bulk crit- 
ical current per unit area a 2 is /bulk _ iQg en A/hkp, 
which is a bit higher than the current density deter- 
mined by the Landau depairing velocity Vd = A/hkp. 
This stems from the possibility of having gapless super- 
conductivity in 3D at superfluid velocities slightly ex- 
ceedingca Vd (note that kp is direction-dependent for a 
cubic lattice at half-filling; we use the average value over 
the Fermi surface kp rj 2.8a, appearing in the trans- 
port formulas, to compare our critical bulk supercurrent 
density to the expressions that assume a spherical Fermi 
surface and a density of particles n = kp/in 2 ). The junc- 
tion properties are studied here in the low-temperature 
limit at T = 0.01 = 0.09T C (the BCS gap is essentially 
temperature independent below 0.6T C ). At this temper- 
ature, the coherence length of the clean normal metal is 
£at = frvp/2irkT ~ 40a. Since we do not consider inelas- 
tic scattering processes, the dephasing length L$ is larger 
than £/y. Therefore, mm L<p) = £jy determines the 
coherence properties of a single quasiparticle wave func- 
tion inside the normal region. 

The inhomogeneous superconductivity problem is 
solved by employing a Nambu-Gor'kov matrix formula- 
tion for the Green function G(rj, r^, iio n ) between two 
lattice sites and Tj at the Matsubara frequency iuj n — 
inT(2n+l), 
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FIG. 3. Charge deviation 8n(z) (from half-filling) in 
SINIS junctions characterized by different Fermi level mis- 
matches AEf = Ep — Ep. The N interlayer consists 
of: (a) 20 normal planes (b) 2 normal planes. The inset 
shows that distributions of Sn(z) for different AEf can be 
rescaled to a single curve after multiplying them by the ra- 
tio (A_B_F) re f I 'AEf , where {AEf)ycI = 1.0 is chosen as the 
reference distribution. 

The corresponding local self-energy is given by the ma- 
trix 

±(T t) iu, n ) = ( ffi^l ti r ;' ZWn \V (3) 

11 ' \ <f>*(Ti, iu n ) -£*(ri,iw n ) J w 

The diagonal and off-diagonal (i.e., normal and anoma- 
lous) Green functions are defined, respectively, as 



G(Ti,Tj,iu3 n ) 



drexp^ n r)(T T c JCT (r)4(0)), (4) 



F(ri,Tj,iw n ) = - J dr exp(iw„r) (T T Cj T (r)cj4,(0)) , (5) 
o 

where T T denotes time-ordering in r and (3 = l/T. The 
self-energies and Green functions are coupled together 
through the Dyson equation, 
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FIG. 4. Electron filling n(z) [dashed line] and charge devi- 
ation Sn(z) [solid line] of a SINIS Josephson junction with 
a thickness L — 20a and the N interlayer chosen to approxi- 
mate a highly doped semiconductor. The charge deviation is 
measured with respect to the equilibrium filling in the bulk, 
ns = f in the superconductor and njv = 0.01 in the normal 
region. The Fermi energy mismatch AEf = Ep — E F be- 
tween the N and the S is: (a) AE F = 1.0, (b) AE F = -1.0, 
(c) AE F = -3.0, and (d) AEf = -10.0. The charge profile is 
virtually independent of temperature, both above and below 
T c . 

G(ri,Tj,iuj n ) = G°(r i ,r j ,iu; n ) 

+ ^ G°(r i ,ri,iuj n )t(ri,iuj n )G(ri,r : j,iu! n ), (6) 
i 

where the local approximation for the self-energy, 
E(rj,rj,iw„) = £( r i> ioj n )5ij, is assumed. In the HFA 
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FIG. 5. Evolution of charge accumulation region in a 
SINIS junction, with a highly doped semiconductor as the 
N interlayer (n N = 0.01, n s = 1.0, AEf = 0.2), upon 
crossing T c of the S by going from T = 0.2 = 1.8T C to 
T = 0.01 = 0.09T C . The inset shows the relative change 
of the charge deviation [5n(z)T=o.2 — Sn(z)T=o.oi]/Sn(z)T=o.2 
for two different charge redistributions: AEf = 0.2 and 
AEf = -0.2. 

for the attractive Hubbard model, the local self-energy 
is found from the local Green function by 

T,(ri,iu) n ) = UiT^2G(r l: r i: iuj n ), (7) 

and 

0(r,-,iw„) = -UiT^2F(ri,Ti,iu n ). (8) 

The self-energy is time- independent because the interac- 
tion is instantaneous and we use the HFA (i.e., retarda- 
tion effects in the superconductor are neglected). The 
noninteracting Green function, G (ri,rj,iu> n ) is diago- 
nal in Nambu space, with an upper diagonal component 
given by 

G°(ri,Tj,iu n ) = / <2 3 k- ■ . (9) 

J lUJ n + [I - £ k 

As all sites within a plane are identical, the self-energy 
need only be calculated once for each of the planes, while 
it is allowed to vary from plane to plane. 

We work with Green functions G a ^/3(iu) n ,k x ,ky) rep- 
resented in a mixed basis, which is defined by the 
two-dimensional momenta (k x ,k y ) and the (discrete) z- 
coordinate of the plane a = zi/a. This follows af- 
ter the initial 3&, problem is converted to a quasi-one- 
dimensional onec2l by performing a Fourier transforma- 
tion within each plane (where the junction is translation- 
ally invariant) and retaining the real-space representation 
for the z-direction of the inhomogeneity. For the local in- 
teraction treated in the HFA, computation of the Green 
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function reduces to inverting an infinite block tridiago- 
nal Hamiltonian matrix in real space. The Green func- 
tions are thereby evaluated as a matrix continued frac- 
tion (technical details are given elsewhereE3o) . The fi- 
nal solution is fully self-consistent in the order parameter 
|A(z)|e^ < - z ' ) inside the part of the junction comprised of 
the N region and the first 30 planes inside the supercon- 
ducting leads on each side of the N interlayer (see Fig. [I]). 
The self-consistent region is long enough because |A(z)| 
heals to its bulk value over just a few coherence lengths 
£o- Our Hamiltonian formulation of the problem and 
its solution by this Green function technique is equiva- 
lent to salving a discretized version of the Bogoliubov-de 
GennesEa (BdG) equations in a fully self-consistent man- 
ner, i.e. by determining the ofltdiagonal pairing potential 
Aj in the BdG Hamiltoniano after each iteration until 
convergence is achieved. The tight-binding description of 
the electronic states also allows us to include an arbitrary 
band structure or unconventional pairing symmetryLJ 

In conjunction with the self-consistent solution of the 
superconducting part of the problem, we have to self- 
consistently solve the electrostatic problem. Although 
both the S and N are half-filled in most of our calcula- 
tions (i.e., there is no mismatch in the Fermi wave vec- 
tor), shifting the bottom of the N band leads to a differ- 
ence in their Fermi levels. This generates a redistribu- 
tion of electrons around the SN interface when these are 
brought into contact. The resulting non-uniform electric 
field can be described by a potential V(z) (for simplicity, 
we use the label z having in mind a discrete z+ coordinate 
at a particular site i) which varies in the transition layer 
around the SN boundary with a thickness of 2c?. In the 
region z > \d\ the following condition is satisfied 



Ep + V(z < 



d) = E$ + V(z >d)= n, 



(10) 



in order to ensure a constant electrochemical potential 
fi throughout the system in equilibrium. The solution 
which satisfies this equation is usually simplified^ to 
V(z) = Vq5(z) + fj, — v(z), where v[z) is a monotonic 
function of z equal to Ep for z < —d or Ep for z > d 
(this allows one to formulate quasiclassical equations in 
the region \z\ > d). Here we treat the contact be- 
tween the S and N in a fully microscopic fashion: start- 
ing from the Hamiltonian (|lj), a Fermi level mismatch 
AEp — Ep — Ep, and assuming a screening length lp> of 
a few lattice spacings, we find the charge redistribution 
around the contact, as well as the corresponding clas- 
sical electrostatic potential generated by them. Thus, 
our technique can treat arbitrary spatial variation of the 
(lattice) Green functions, superconducting order param- 
eter, and electrostatic potential. This includes the region 
\z\ < d, where we find a sharp increase of V(z) but never 
as sharp as the (unphysical) delta function. 

Since our SINIS multilayer structure is translation- 
ally invariant in the transverse direction, each infinite 
plane has a uniform surface charge distribution Sn(z)a 
which generates a homogeneous electric field E(z) = 



8n{z)a/2eae oa pointing along the z direction (eoo is the 
relative dielectric constant of the ionic lattice). The 
quantum-mechanical part of the electrostatic problem en- 
tails determining the local electron density n(zi) = n(z) 
(filling) at each site of a given plane a — z/a 

oo 

n(z) = k B T^2 / p 2D {e X y)\va.G aa {iLo n ,£ xy )de X y, 

(11) 

where e xy = — 2f [cos(k x a) + cos(fcj,a)] is the in-plane ki- 
netic energy for the transverse momentum (k x ,k v ), and 
p 2D (e xy ) is the two-dimensional tight-binding density of 
states on a square lattice (which is used for the sum over 
momenta parallel to the planes). The corresponding elec- 
tric potential is determined classically from the "charge 
deviation" Sn(z) — n{z) — n (n is the average filling in 
the bulk, tin or rig) 



5V{z) = - 



eaSn(z')\z — z'\ 
2eo £ oo 



(12) 



This must be summed over all planes to give the on-site 
potential V(z). Therefore, the small induced charge im- 
balance 5n(z) — N(fi)e6V(z) satisfies (in a corresponding 
continuous system) 



d e 2 aN(n) 

— dn(z) = on(z) 

dz 2e eoo 



(13) 



where N{fi) is the total density of states at the chemical 
potential fi. This is integrated to give the distribution of 
the screened charge 



Sn(z) = Sn(zo) exp 



2e e c 



(z - z ) 



(14) 



which decays exponentially on a length scale set by the 
Debye screening length 



Id = 



e 2 aN{n) 



2e n e 



£ oo 



(15) 



Thus, the screening length is determined by N(/j,) and 
too (for example,U £oo = 20 — 30 and Id = 5 — 10 A in 



high-T c superconductors). We choose ej; 
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which leads to Id ~ 3a. The self-consistency in the elec- 
trostatic problem is required because V(z) enters into the 
computation of the Green function as a diagonal poten- 
tial in the Hamiltonian ([[]). The solution has converged 
when the potential is consistent with the charge distribu- 
tion (pj]) determined from the Green function. Although 
this seems like a cumbersome computational task, the 
potential around the SN boundary barely changes when 
equilibrium Josephson current flows. Thus, the electro- 
static part of the problem converges rapidly since the 
potential found in the solution at one phase gradient is 
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a good initial guess for the iteration scheme at the next 
superconducting phase gradient. 

The density of electrons n(z) on each site in a given 
plane (at zero supercurrent) is plotted as a function of 
the junction thickness for AEp = 3.0 in Fig. [| The 
charge deviation Sn(z) from half-filling n$ = n^ = 1 and 
the corresponding electrostatic potentials are (approxi- 
mately) symmetric around the SN boundary for thick 
enough junctions, as shown in Fig. || Strictly speak- 
ing, only such symmetric distributions should be denoted 
"screened dipole layers" . For thinner junctions, where 
the screening of the excess charge does not heal n(z) to 
its equilibrium value, charge is depleted from the N in- 
terlayer. The example of this behavior is the L = 2a 
junction in the lower panel of Fig. |[ It leads to a non- 
monotonic resistance as a function of junction thickness 
L at fixed AEp (see Sec. [V). Thus, the charge effects be- 
come essential in short-coherence length superconducting 
junctions with thicknesses L < 2L«, which are encoun- 
tered in high-T c grain boundaries E3 

The interesting feature of the 8n(z) profiles for the half- 
filled S and N is that they can be approximately rescaled 
to a single reference distribution [set by (AEp) re f\ after 
multiplying each of them by the ratio (AEp) Ie {/ AEp, 
as shown in the insets of Fig. ||. We believe this oc- 
curs because the noninteracting cubic density of states 
is nearly constant close to half filling. The scaling be- 
comes essential in computing the properties of junctions 
with large AEp since one can use the scaled potential 
profile computed at a smaller AEp as the initial guess 
in the iteration procedure. Since the n$ — rijv = 1 
case has a higher degree of symmetry, we also per- 
form calculations for = 0.01 (which approximates 
a doped semiconductor) in the normal region and half- 
filling ns = 1.0 in the superconductor. The result is 
shown in Fig. ||. Here the scaling of the 8n{z) distribu- 
tion does not work as well (because the density of states 
has strong variation with energy). In addition, we find 
that the charge deviation is nonsymmetric, and yields a 
different 5n(z) for positive and negative AEp [for sym- 
metric filling ns = tin the two distributions are simply 
related as Sn(z)\^AE F — —Sn(z)\AE F ]- We also investi- 
gate the temperature dependence of the distributions of 
uncompensated charge for T — 0.2 (the chemical poten- 
tial in the bulk N is fi = —5.566 for njv = 0.01) and at 
T = 0.09 (which is close to T c = 0.11). In both cases 
ns = n N = 1-0 and ns ^= tin = 0.01 we find 5n(z) to be 
practically temperature independent (e.g., the change is 
at most 5% around the SN boundary) for AEp shown in 
the previous figure. This feature is exploited in Sec. IV 
to calculate the normal state resistance of our junctions 
from an imaginary axis computation of the charge and 
potential profile in the superconducting state. However, 
for ns 7^ n^ — 0.01 and small \AEp\ ~ 0.2 a large 
change in the magnitude of Sn(z) is observed when go- 
ing from T = 0.2 > T c to T = 0.01 < T c , as shown in 
Fig. |^. Similar phenomenon has been found in the re- 
cent Raman studies^ which show a substantial change 



in the thickness of the charge accumulation layer at the 
interface between Nb and InAs, as Nb undergoes a su- 
perconducting transition and proximity effects develop 
in the InAs layer. This would point to a proximity ef- 
fect influenced screening length, which cannot be seen 
in our local (Thomas-Fermi) screening theory containing 
only two parameters which determine Id'- £oo, which is 
fixed in our calculations, and the density of states N(p) 
which can be modified by the proximity effect. Our ob- 
servation of the change in the charge concentration above 
(T ~ AEp) and below T c , without a palpable change in 
the screening properties, suggest that effects beyond the 
simple screening theory (e.g., nonlocal-|Screening which 
becomes important in low filling casesa) probably have 
to be taken into account to understand this experiment. 



III. SELF-CONSISTENT EQUILIBRIUM 
PROPERTIES OF SINIS JUNCTIONS 

We first provide an insight into the microscopic prop- 
erties of these junctions which are determined by the 
proximity effect that affects the critical current (in non- 
self-consistent calculations such effects are taken into ac- 
count only through spme effective phenomenological sup- 
pression parametem). They are encoded in the self- 
consistently computed variation of the amplitude and 
phase of the order parameter A(z) = |A(0)|e^^ z - ) in the 
S or pair amplitude F(z) = |F(z)|e^W in the N. These 
are related to each other inside the S by 



A(z) = -U(z)F(z), 



(16) 



where F(z) is obtained as the equal-time limit of the local 
anomalous Green function introduced in Sec. |F] 



F(z) =F(z hZi 



H 



(17) 



Although 4>{z) and F(z) are not directly measurable, 
they are important for understanding the superconduc- 
tivity in inhomogeneous structures. Examples include 
the proximity effect in the N and the depression of |A(z)| 
(compared to its bulk value) on the S side of a SN bound- 
ary ( "inverse proximity effect" ) . Since the critical current 
of the junction is determined by A(z) at the SN bound- 
ary, the study of F{z) throughout the junction gives di- 
rect insight into how self-consistency affects the transport 
properties (analytical approaches usually assume a step 
function for A(z)|, which is applicable only for a limited 
range of junction parameterst3) . The non-zero value of 
F{z) inside the superconductor results from the attrac- 
tive pairing interaction U(z) ^ [which also gives rise 
to the non-zero order parameter A(z)]. In the normal 
metal, U(z) = and the gap vanishes, but F (z) can be 
non-zero due to the proximity effect. Therefore, it is more 
meaningful to plot F(z), which is a continuous function 
throughout the junction. Inside the S, F(z) should be 
understood as just A(z)/[—U(z)]. The superconducting 
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FIG. 6. Pair amplitude F(z) at zero supercurrent flow 
in SINIS junctions of different thicknesses L = TVjva. 
The double-barrier structure arises from the charge inhomo- 
geneity (see Fig. ^ induced by the difference between the 
Fermi energies of the normal metal and superconductor: (a) 



AE f 



E§ 



1.0, and (b) AEf = E£ 



E F = 3.0. 



correlations are imparted to the N region which is in 
contact with the S'^-They are described quantitatively by 
the pair amplitude^ F(z) dl7|). Because of the transla- 
tional symmetry of the junction in the transverse direc- 
tion, F(z) is constant within the plane, and changes from 
plane to plane along z-axis. The scale over which F[z) 
changes exponentially from the SN interface to zero in 
the bulk of the N is set by the normal metal coherence 
length £jy. However, as T — > the length ^ diverges 
and the exponential decay of F(z) crosses over to a slower 
power-law decay £like 1/z at T = 0, inside a N described 
by a Fermi liquico). 

Although this description of the proximity effect has 
been used since the eatly days of inhomogeneous su- 
perconductivity studies,c2l it is only recently that meso- 
scopic superconductivity!] has established an explicit con- 
nection to a real-space picture of pairing correlations, 
provided by the— phenomenon of (phase-coherent) An- 
dreev reflect ion That F(z) is non-zero in a normal 
region is equivalent to saying that the electron and an 



Andreev reflected phase-conjugated hole maintain their 
single-particle phase coherence inside the N. Technically, 
this interpretation follows directly from the expression for 
F(z) in terms of-auasiparticle wave functions entering the 
BdG equations.Q In other words, near the SN boundary, 
Andreev reflection mixes electron-like and hole- like quasi- 
particles in the same proportion in which they are mixed 
in the S (where Bogoliubov quasiparticles are a mixture 
of electron-like and hole-like states with weights deter- 
mined by the self-consistency condition) due to purely 
kinematic effects, since the interaction is absent in the 
N. The definition of F{z) from Eqs. (§) and @ in 
the second-quantized formalism, shows that such corre- 
lations can be interpreted alternatively as a condensate 
wave function leaking into the normaLmetal through the 
presence of evanescent Cooper pairs. c3 In the case of a 
Josephson junction, the overlap of two condensate wave 
functions provides a weak coupling between the super- 
conducting leads, while insuring the global phase coher- 
ence and equilibrium current flow (i.e., the DC Joseph- 
son effect) for the time-independent phase difference be- 
tween them. Thus, the two apparently different pic- 
tures of the Josephson effect in weak links (leakage of 
Cooper pairs versus Andreev reflection induced transfer 
of Cooper pairs) are in fact two facets of the same phe- 
nomenon. 

We first show two examples of F(z) computed self- 
consistently for vanishing supercurrent throughout the 
SINIS junction with ns — tin — 1. Figure ^| plots the 
scaling of F{z) with the junction thickness for the I layers 
at the SN boundary being SDLs whose height is deter- 
mined by AEf = 1.0 or AE F = 3.0. The shape of F(z) 
evolves with the thickness, as well as with the height of 
the double-barrier. This second point is demonstrated in 
Fig. ^ where we fix L and vary the strength of the SDL 
barrier. Here one would expect the evolution of F(z) to- 
ward a step function, which then justifies the use of rigid 
boundary conditions for strong enough scattering at the 
/ barriers.tj However, we find a non-monotonic change 
in the shape of F(z): the influence of a SDL on the order 
parameter A(z) is first reduced with increasing AEp, 
but then leads to a depressed A(z) near the boundary 
for a strong charge imbalance generated by AEp = 10.0. 
Since our previous results for a SINIS junction having a 
strong on-site Coulomb potential, confined within a sin- 
gle plane, exhibit a step function UkeEa A(z), the effects 
observed here can be attributed to the finite spatial ex- 
tent of the SDL. Moreover, we find that the step function 
(up to tiny oscillations near the boundary) for A(z) does 
develop in the special case of low filling in the N region, 
like n^v = 0.01, and a small mismatch |A£j?| < 1. A 
specific example of this behavior (compared to the case 
with the same parameters, but with a negative AEp) is 
shown in Fig. 

In the short junction case, the oscillations of A on the 
scale of Xp are observed for large enough AEp. In this 
case, as discussed in the previous section, the junction is 
too thin for the distribution of charge to heal to its 
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FIG. 7. Pair amplitude F(z) at zero supercurrent flow in 
SINIS junctions of thickness L — 20a characterized by differ- 
ent heights of the SDL barriers. The double-barrier structure 
(charge inhomogeneity from Fig. ^) arises from the difference 
in Fermi energies of the normal metal and superconductor, 



AE F = Ep 



equilibrium value. The charge depletion inside the N 
brings it close to an insulating state. While oscillations 
on the scale of A_f have been observecO in similar self- 
consistent calculations at T = Q (and attributed to the 
mesoscopic coherence of a single particle wave function), 
here it appears that they are a property of the super- 
conducting interface which terminates at an "insulator" 
(this is also exhibited by a thick junction with small 
in Fig. ||). We have recently found such behavior, in its 
most pronounced form, in the case of SIS junctions, with 
/ being a correlated insulator e3 In the self-consistent 
calculations one can also observe how F(z) evolves, be- 
coming smaller inside the N region, as the phase change 
across the interlayer is increased and the Josephson cur- 
rent approaches I c . An example of such an effect due to 
self-consistency is shown in Fig. ||. 

When self-consistency is satisfied, the phase of the or- 
der parameter \A(z)\e l ^ z ^ is not a constant inside the S 
leads (see also Sec. |V|) because a phase gradient is needed 
to support the current in the S ensuring current conserva- 
tion throughout the junction. Thus, the change of phase 
from plane to plane has to be extracted from the self- 
consistent solution for \F(z)\e l ^ z \ It can be expressed 
as the sum of a linear term and a "phase deviation" term 

s<Kz) 



dz 



■6<t>(z), 



(18) 



bulk 



where the distance z is measured from the origin along 
the z-axis. The linear term is determined by the phase 
gradient (dip / 'dz)buik which is set as the boundary condi- 
tion in the bulk of the superconductor. The non-trivial 
information contained in <j)(z) is revealed by plotting 
8<j)(z). The overall phase <fi(z) increases smoothly and 
monotonically across the self-consistent region. We plot 
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FIG. 8. Pair amplitude F(z) at zero supercurrent flow in 
SINIS junctions with thickness L — 20a. The double-barrier 
structure arises from the inhomogeneous charge redistribu- 
tions, plotted in Fig. ^, induced by AEf = Ep — Ep between 
the normal region with njv = 0.01 and the superconductor at 
half-filling ns = 1. 

8<f)(z) for two different junction thicknesses with 
AEf = 1 hi Fig- f[o| In general, oscillations of 5<p(z) 
on the scale of Af are found for moderate AEp and long 
enough junctions (L > 6a). Oscillations, of both the 
phase and the order parameter, were found inside a long 
mesoscopic constriction in previous self-consistent calcu- 
lationsE-3 at zero temperature (that gradually disappear 
with increasing T). Here we see the oscillations of 5<p(z) 
at low temperature (but still L < £at), while the corre- 
sponding |F(z)| (Fig. ||) does not oscillate. 



IV. CRITICAL CURRENTS AND 
CHARACTERISTIC VOLTAGES 

In the self-consistent treatment, equilibrium supercur- 
rent flows through the junction when the phase gradient 
(d(j) / dz)huXk exists in the bulk of the superconductor and 
a total phase change cf> is established across the normal 
region. Therefore, we first find the solution for the bulk 
superconductor in both the absence of a supercurrent and 
in the presence of a supercurrent generated by a uniform 
variation in the order-parameter phase. The uniform 
bulk solution is then employed to provide the "boundary 
conditions" for the junction beyond the region where we 
determine properties self-consistently. Thus, our method 
does not require any assumptions about the boundary 
conditions at the interface between the barrier and the 
superconductor_which follow from the requirements of 
self-consistency.E£l We use current conservation as a strin- 
gent test of the achieved self-consistency in the solution 
for the Green function. Namely, the self-consistently de- 
termined A(z) ensures that Andreev reflection at each 
SN boundary generates supercurrent flow in the S leads 
(besides being responsible for the proximity effect in 
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FIG. 9. Pair amplitude \F(z)\ at different supercurrent 
flows through the SINIS junction as a response to a non-zero 
phase gradient in the bulk and corresponding phase change (f> 
across the N region. The thickness of the junction is L = 20a, 
and the two SDLs are generated by the Fermi energy mis- 
match AEf = 1.0 between the N and S, which are both at 
half filling. Panel (a) is a blow up of the dashed square region 
in panel (b). 

the N discussed in Sec. O. Thus, the fulfillment of the 
self-consistency condition (|16[) means that the "source 
term" (on the right-hand side) vanishes in the equation 
of motion for the charge density operator hi 

dhi ^-^ T 2ie 



at 



•2> 



(A^c^-A?^)), (19) 



thereby recovering current continuity at every site (Zy is 
the current between two neighboring sites). When the 
current inside the superconductors is small, e.g., due to 
the geometrical dilution of a weak link with a junction 
area much smaller than £q, or when the junction rosis» 
tance is dominated by a large interlayer resistanceJ§E2l 
one usually neglects the supercurrent flow and corre- 
sponding phase gradient in the bulk superconductor nec- 
essary to support it. Strictly speaking, such approaches 
violate current conservation.t3L3 Inasmuch as our 5* and 
TV layers have the same area, I c //,!? ulk can be close to one 
for thin junctions with weak SDLs at small AEp. In such 
cases, current flow affects appreciably the superconduct- 
ing order parameter [i.e., F(z) both inside and outside 
of the N, cf. Sec. till and a self-consistent treatment be- 



comes necessary (as is the-case for the critical current of 
the bulk superconductora) . Because of the presence of 
a phase gradient inside the S, the simple picturetl of an 
equilibrium current being related to the phase difference 
4>l — 4>r between the left and right S leads (where 4>l 
and <f>R are constant within the leads) is not applicable. 
Nevertheless, the solution for the current turns out 
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FIG. 10. Phase deviation 5(j>{z) [the total phase change 
along z-axis is the sum of the bulk phase gradient and 5<f>(z), 
Eq. (E3)], within the self-consistently modeled part of the 
SINIS junction, at different supercurrent flows (the last 
curve corresponds to I c ). The thicknesses of the junctions 
are: fa) L = 2a, and (b) L = 20a, and the double-SDL-barrier 
(Fig. ph is induced by the Fermi energy mismatch AEf = 1.0 
between the iV and S. 

to be uniquely parameterized by a single quantity 
which j-can be taken as the phase change access the N 
regionc3 (the other option is the phase offsefO which is 
related to the phase change by a nontrivial scale transfor- 
mation). In a discrete model like ours, a convention has 
to be introduced on how this change is extracted from 
<f>(z) in Eq. ([l8|). The thickness of the junction is defined 
to be the distance measured from the point Zl, in the 
middle of the last S plane on the left (at zf) and the 
first adjacent N plane (at = zf + 1), to the middle 
point zr between the last iV and first S plane on the 
right (cf. Fig. |l|). Since (f>(z) is defined within the planes, 
we set 4>{zl) = [^(^f) + 4>( z l)]/^ to be the phase at zl, 
and equivalently for cf>(zn). The phase change across the 
barrier is then given by 



C,)-o(, L ) = I(g 



bulk 



(z L ), 
(20) 
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FIG. 11. Scaling of the current- phase relation I(<j))/I c with 
the thickness of the SINIS junctions where the SDLs are 
determined by: (a) AE F = 1.0, and (b) AE F = 3.0. Note 
that the phase change <f> c at the critical current I c = I(<j} c ) 
varies nonmonotonically with the junction thickness, as shown 
in Fig. [12J The standard I((f>)/I c = sin</> dependence in the 
SIS tunnel junctional is plotted as a reference only (which is 
analytically predicted for SINIS junctionsrwith small barrier 
transparency at high enough temperatures"). 

for a junction of thickness L. The current-phase re- 
lation I((f>) is obtained by computing the current for 
a fixed bulk phase gradient and associating this value 
with the phase change across the N region, which is ex- 
tracted from the self-consistent pair amplitude F(z) = 
\F (z)\e l ^ z \ On the lattice, transport is described by 
the current across a link between two adjacent planes a 
and a + 1. This current (per a 2 ) is obtained from Ihc 
Green function connecting two neighboring planes asE3 



oo 

2eat f 2n 

iq,a+i = —^~k B T ) j J p (e xy ) 



'dm [G a ,a+l{iu n , £xy)] d,£ xy . 



(21) 



The first iteration in our self-consistent algorithm usually 
gives a current which is smaller inside the N than in 
the S region. The iteration cycle is completed when the 
current is constant throughout the junction. The only 
approximation invoked here is the presence of a (typically 
small) discontinuity in the supercurrent at the 



FIG. 12. Phase change C at the critical current I(<j)c) ~ Ic 
plotted against the thickness of the SINIS junctions, where 
the double-SDL-barrier is determined by AEf = 1.0 (circles), 
AEf = 3.0 (squares), and AEf = 5.0 (diamonds). 

bulk-superconductor/self-consistent-superconductor 
interface. We find that the superconducting order has 
always healed to its bulk value at that point. However, 
sometimes there can be a jump in cj)(z) at this boundary 
when one nears the critical current. This discontinuity 
in the phase corresponds to a breakdown of current con- 
servation at the this interface (it can become large for 
large AEp and a thick junction, especially when one lies 
on the decreasing current side of the current-phase dia- 
gram). The critical current I c of the junction is reached 
when the planes with the lowest pair amplitude |-F(#)| 
(which are located in the center of the N) can no longer- 
support the necessary phase gradient to maintain current 
continuity. 

The scaling of the shape of the current-phase relation 
with the junction thickness is plotted in Fig. |ll| for dif- 
ferent AEp- We find large deviations from the usual 
sinusoidal I(cf>) dependence for thin junctions and mod- 
erate heights of the SDL barriers. While in such feafPj 
(and at low temperatures) analytical predictionsoliElJ 
also give non-sinusoidal I(cf>) , our "critical" phase change 
(j) c [I{4>c) — Ic] is always below the analytical prediction 
4> c w 1.86 (Fig. naL-which can be attributed to the effects 
of self-consistencyEil (the other important distinction is 
that SDLs are spatially extended barriers). For thicker 
junctions, with high SDL barriers (and at high enough 
temperatures) the recovery of the usual SIS junction 
I(<j)) — I c sin(j) current-phase relation is predictedB Here 
we find a current-phase relation /(</>) which is close to si- 
nusoidal in the thick junction limit [(a) panel in Fig. |ll| ], 
or in thin junctions with high SDL barriers [(b) panel in 
Fig. ^l| . The corresponding critical current densities as a 
function of junction thickness are plotted in Fig. [l]. For 
large AEp, I C {L) is non-monotonic because of the special 
role played by the barriers formed in the junctions with 
L < 21 p>. When SDLs are completely screened inside the 
thick N interlayers, the decay of current is determined 
just by the exponential decay of the proximity 
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FIG. 13. Semilogarithmic plot of the critical current (per 
unit area a 2 ) as a function of junction thickness (i.e., number 
of planes iVjv =1-60 inside the N interlayer) of the SIN IS 
junction for different SDLs determined by AEf = Ep — Ep. 
Both the S and N are at half-filling in the bulk. 

coupling between the S leads through the clean nor- 
mal interlayer (the resistance of these junctions is also 



practically independent of L, see Fig. 14). For exam- 

P le >r4 
tingP 
in Fig 



. characteristic decay length, extracted from fit- 
I ALPexp[-L/£' N } to I C (L) for AE F = 1.0 case 
13], is £' N « 35a ~ £n (with p = 0.3). For larger 



AEf, and long enough junctions to ensure monotonic 
decay of I C {L), £' N appears to be shorter. 

We use a Kubo linear response formalism to determine 
the normal state resistance Rjy = 1/Gm- Kubo theory is 
formulated in terms of the non-local conductivity tensor 



j(r,f) 



ydr'a(r,r»-E(r'^), 



(22) 



which relates the current density j(r) to the electric field 
E(r) through a non-local Ohm's law (at finite frequency 
v these are the respective Fourier components) . Its phys- 
ical meaning is obvious — it gives the current response at 
r due to an electric field at r'. Although an external 
electric field induces charges (and corresponding poten- 
tials) to linear order, the linear transport properties, like 
cr(r,r'), are found as the response to an external field 
only. This is because the current response to this inho- 
mogeneous field-Xexternal + induced) is already beyond 
linear responseo Thus, only equilibrium screening has 
to be included in the Hamiltonian used— to compute the 
Green function entering a(r, r') below. EHl This makes it 
possible to use the potential generated by the charge dis- 
tributions 8n(z) (discussed in Sec. ||), which is computed 
from the imaginary axis calculations, as an on-site fixed 
potential in the equilibrium Hamiltonian ([lj) . In this way 
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FIG. 14. Semilogarithmic plot of the normal state resis- 
tance (multiplied by the unit area a 2 ), as a function of the 
SIN IS junction thickness, for different SDLs determined by 
AE F = Ep — E F . Both the S and N are at half-filling 
in the bulk. Note that junctions with thickness smaller 
than 21 d have charge depleted for large enough AEf causing 
the resistance to change non-monotonically as a function of 
L. The Sharvin point contact resistance of the clean SNS 
junction (corresponding to AEf = in our structure) is 
R sh a 2 = l.58ha 2 /2e 2 

the potential of the SDLs enters the resistance calcu- 
lation through the Green functions in Eq. computed 
by real-axis analytic continuation. The DC conductance 
of the sample of volume is expressed through er(r, r') 

as 

^ = 7)2/ drE(r)-j(r) 

n 

= -L J dr dr' E(r) • a(r, r') • E(r'), (23) 
n 

where E(r) is the local field inside the sample and V is the 
externally applied voltage. Because of current conserva- 
tion requirements on the form oft3 a(r,r') ! _it is possible 
to use arbitrary electric field factors in Eq. ( |23j ) [including 
a homogeneous field E = V/L]. 

Since our system is effectively one-dimensional (in real 
space) we need to calculate the longitudinal component 
in the z-direction (perpendicular to the uniform planes) 
of <x(zi, z[). In a lattice model like ours, the relevant com- 
ponent of this tensor, a zz {a, (3), is given by (neglecting 
vertex corrections) 



5 (a,/3) 



-1 (eatf 



P xy)ds X y 



2^ 
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x [lmG at/3+1 (uj,e X y)ImG g ^ +1 (uj,e X y) - 
Im G a ,p(u, e xy ) Im G a+ha+1 (uj, e xy )] 
x [cosh 2 (cj/2k B T)}- 1 , (24) 

where f{u>) is the Fermi-Dirac distribution function. We 
first find the self-consistent solutions for the system in 
the normal state, with no current flowing, by setting the 
order parameter to zero on all planes. These solutions are 
then employed to calculate the Kubo tensor (24). The 
self-energy of the planes outside the interlayer contains 
only a constant real part, as the calculation is carried out 
within the HFA. Given the set of local self-energies, the 
Green functions which couple any two planes are readily 
found, for any momentum parallel to the planes. 

The conductance Gm (per unit area a 2 ) of the lattice 
system is obtained from the discretized version of (|22 



—5- = [R N a ) 



a,f3 



(25) 



as the sum of the components of the non-local Kubo 
conductivity tensor ..— . Thus, although one can find the 
inhomogeneous fielcES Eg « +1 (across all links connect- 
ing planes /3 and (3+1) by inverting the discretized 
version of Eq. (||), I a , a+1 = aj^g Czz(a, 0)E b ,b+\ for 
Ia,a+i a constant throughout the system, the final ex- 
pression for the conductance does not contain this field. 
The normal state resistances calculated in this frame- 
work are plotted in Fig. [l4|. In thin junctions and for 
large enough AEp, a charge depletion layer arises inside 
the TV which leads to a non-monotonic behavior of 
(e.g., _Rjv increases sharply for L = 2a and AEp = 3.0, 
or L — 3a and AEp = 5.0). On the other hand, 
for small enough AEp < 1.0 the conductance is only 
slightly changed from the Sharvin point contact conduc- 
tancald of a ballistic SNS junction per unit area a 2 , 



R N a 2 = [{2e 2 /h){k 2 



1.58ha 2 /2e 2 . Therefore, 



comparison of Fig. [13] and Fig. |lj shows that the SDL 
depresses the current substantially, while only weakly in- 
creasing the resistance. This reduces the I c Rn product, 
plotted in Fig. [l^, thus showing that charge accumulation 
layers are detrimental to junction performance in elec- 
tronics circuits. This is further confirmed by the fact that 
I c Rn m most of these junction is below the product of 
the bulk critical current and the Sharvin point contact re- 
sistance /,!? ulk i?sh = 1.45A/e, which is the upper limit of 
the characteristic voltage in a clean SNS weak link (the 
SNS junction made of the same S leads as studied here, 
but with a dirty N interlayer ^xhibits I c Rn > 2k ulk i? sh , 
for some range of parametersEj) . Therefore, the SDL in- 
duced scattering on a SN boundary is one of the mech- 
anisms which can account for the low I c Rn products 
observed in experimentsa on nominally ballistic short 
SNS junctions (where Rn, being determined by the thin 
charge layer only, does not scale with L just like what 
happens in ballistic conductors). One way to test this 
conjecture is to use electron holography to map out 
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FIG. 15. Product of the critical current 7 C and the nor- 
mal state resistance -Rjv as a function of the SIN IS junction 
thickness, for different SDLs determined by AEf ~ Ep —Ep. 
Both the S and N are at half-filling in the bulk. The I c Rn 
is always below the product of the bulk critical current and 
Sharvin point contact resistance 7^ ulk ^Sh = 1.45A/e (dashed 
line), except in the case L — 3a(~ Id) with AEp = 5.0. 

the charge profile near the SN interface of these ballis- 
tic junctions. 



V. CONCLUSIONS 

We have studied the influence of a charge imbalance, 
that arises at the boundary between a short coherence 
length superconductor and a normal metal (due to Fermi 
energy mismatch) on the equilibrium properties of a 
SINIS Josephson junction (where the S and N layers 
are of the same width). The screening length is large 
enough to generate a spatially extended charge redis- 
tribution that allows us to examine the interplay be- 
tween the charge layer formation and superconductivity 
(characterized by the coherence length comparable to the 
screening length) near the SN boundary and in the 2V" 
interlayer. This resembles the charge redistribution on 
the grain boundaries of a high-T c superconductor. At 
half-filling in both the S and 2V, the charge distribution 
and its potential are symmetric (screened dipole layer), 
and can be rescaled to a single one determined by some 
reference Fermi level mismatch AEp . When charge con- 
centration in the N is a hundred times smaller than in 
the 5, we find a proximity effect induced change in the 
charge redistribution generated by a small Fermi level 
mismatch upon moving from T > T c (where T is of the 
order of AEp) to T < T c . 

The step-function-like order parameter (which is used 
in non-self-consistent approaches) is recovered only in the 
case of a low charge density in the N (compared to the 
filling in the S) and a small mismatch \Ep — Ef.\ < 1. 
The SINIS junction exhibits unusual properties when 
its thickness is comparable to the screening length. While 
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the charge layer leads to a depression of the order param- 
eter near the SN boundary, and thereby the junction 
critical current, it influences the normal state resistance 
in a much weaker fashion. Therefore, the I c Rn product, 
relevant for digital electronics application, is reduced. 
This points out that such space-charge layers should be 
avoided to optimize junction performance ancLincrease 
the critical current in high-T c superconductors^ 
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